void Write()

//Write histos with all b coefficients into the file Hists.root

{
    FileData = new TFile("MergedOutput.root");
    FileMC =  new TFile("OutputPythia.root");
    FileMCPhojet = new TFile("OutputPhojet.root");

    //TString ListNames[15] = {"StandardCuts","Vertex3","Vertex5","Vertex10","TPC40","TPC100","ITS0","ITS4","DCA01","DCA20","ITSTPCRefitFalse","TPConly","StandardITSTPC","StandardITSTPCfalse","StandardITSTPCCrossedRows"};

    TString ListNames[13] = {"StandardCuts","Vertex5","Vertex10","TPC40","TPC100","ITS0","ITS4","DCA01","DCA20","ITSTPCRefitFalse","TPConly","StandardITSTPC","StandardITSTPCCrossedRows"};

    TH2D *HistBData = new TH2D("fHistBData","2D hist for all bData corr coefficients", 16, 0, 16, 15, 0, 15);
    TH2D *HistBMC = new TH2D("fHistBMC","2D hist for all bMC corr coefficients", 16, 0, 16, 16, 0, 16);
    TH1D *HistBMCtrue = new TH1D("fHistBMCtrue","hist for bMCtrue corr coefficients", 16, 0, 16);
    TH1D *HistBMCtruePhojet = new TH1D("fHistBMCtruePhojet","hist for bMCtrue corr coefficients Phojet", 16, 0, 16);
    TH2D *HistB = new TH2D("fHistB","2D hist for all b corr coefficients", 16, 0, 16, 16, 0, 16);

    for(Int_t i = 1; i <= 13; i++)
    {
        HistBData->GetYaxis()->SetBinLabel(i,ListNames[i-1].Data());
        HistBMC->GetYaxis()->SetBinLabel(i,ListNames[i-1].Data());
        HistB->GetYaxis()->SetBinLabel(i,ListNames[i-1].Data());
    }

    HistBData->GetYaxis()->SetBinLabel(14,"Correlator");
    HistBMC->GetYaxis()->SetBinLabel(14,"Correlator");
    HistB->GetYaxis()->SetBinLabel(14,"Correlator");

    HistBData->GetYaxis()->SetBinLabel(15,"FB-BF symmetry");
    HistBMC->GetYaxis()->SetBinLabel(15,"FB-BF symmetry");
    HistB->GetYaxis()->SetBinLabel(15,"FB-BF symmetry");

    HistB->GetYaxis()->SetBinLabel(16,"Phojet correction");
    HistBMC->GetYaxis()->SetBinLabel(16,"Phojet");

    HistBData->GetXaxis()->SetTitle("#eta configuration");
    HistBMC->GetXaxis()->SetTitle("#eta configuration");
    HistBMCtrue->GetXaxis()->SetTitle("#eta configuration");
    HistBMCtruePhojet->GetXaxis()->SetTitle("#eta configuration");
    HistB->GetXaxis()->SetTitle("#eta configuration");


    TList *ListMCtrue = FileMC->Get("MCtrue");
    TList *ListMCtruePhojet = FileMCPhojet->Get("MCtrue");

    for (Int_t i=1;i<=16;i++)  //loop over different Eta configurations
    {
        THnSparseF* SparseMCtrue= (THnSparseF*) ListMCtrue->FindObject(Form("fPhi2piEta%d",i));
        TH2D* CorrPlotMCtrue=SparseMCtrue->Projection(6,7,"");
        TH1D* CorrFuncMCtrue=(TH1D*)CorrPlotMCtrue->ProfileX("");
        CorrFuncMCtrue->Fit("pol1","");
        Double_t bMCtrue=CorrFuncMCtrue->GetFunction("pol1")->GetParameter(1);
        HistBMCtrue->SetBinContent(i,bMCtrue);

    }

    for (Int_t i=1;i<=16;i++)  //loop over different Eta configurations
    {
        THnSparseF* SparseMCtruePhojet= (THnSparseF*) ListMCtruePhojet->FindObject(Form("fPhi2piEta%d",i));
        TH2D* CorrPlotMCtruePhojet=SparseMCtruePhojet->Projection(6,7,"");
        TH1D* CorrFuncMCtruePhojet=(TH1D*)CorrPlotMCtruePhojet->ProfileX("");
        CorrFuncMCtruePhojet->Fit("pol1","");
        Double_t bMCtruePhojet=CorrFuncMCtruePhojet->GetFunction("pol1")->GetParameter(1);
        HistBMCtruePhojet->SetBinContent(i,bMCtruePhojet);
    }

    for (Int_t j=1; j<=13;j++) //loop over different cuts
    {

        TList *ListData = FileData->Get(ListNames[j-1]);
        TList *ListMC = FileMC->Get(ListNames[j-1]);


        for (Int_t i=1;i<=16;i++)  //loop over different Eta configurations
        {
            THnSparseF* SparseData= (THnSparseF*) ListData->FindObject(Form("fPhi2piEta%d",i));
            TH2D* CorrPlotData=SparseData->Projection(6,7,"");
            TH1D* CorrFuncData=(TH1D*)CorrPlotData->ProfileX("");
            CorrFuncData->Fit("pol1","");
            Double_t bData=CorrFuncData->GetFunction("pol1")->GetParameter(1);
            HistBData->SetBinContent(i,j, bData);

        }

        for (Int_t i=1;i<=16;i++)  //loop over different Eta configurations
        {
            THnSparseF* SparseMC= (THnSparseF*) ListMC->FindObject(Form("fPhi2piEta%d",i));
            TH2D* CorrPlotMC=SparseMC->Projection(6,7,"");
            TH1D* CorrFuncMC=(TH1D*)CorrPlotMC->ProfileX("");
            CorrFuncMC->Fit("pol1","");
            Double_t bMC=CorrFuncMC->GetFunction("pol1")->GetParameter(1);
            HistBMC->SetBinContent(i,j, bMC);

        }

        for (Int_t i=1;i<=16;i++)  //loop over different Eta configurations
        {
            if(HistBMC->GetBinContent(i,j))
            {
                HistB->SetBinContent(i,j, (HistBData->GetBinContent(i,j))*(HistBMCtrue->GetBinContent(i))/(HistBMC->GetBinContent(i,j)) );
            }
            else {cout<<"ERROR! bMC = 0";}
        }

    }

    //correlator__________________________________________________________________________________
    for (Int_t i=1;i<=16;i++)
    {
        TList *ListData1 = FileData->Get(ListNames[0]);
        THnSparseF* SparseData = (THnSparseF*) ListData1->FindObject(Form("fPhi2piEta%d",i));
        TH1D *QfQbData = SparseData->Projection(8,"E");
        TH1D *QfData = SparseData->Projection(7,"E");
        TH1D *QbData = SparseData->Projection(6,"E");
        TH1D *QfQfData = SparseData->Projection(9,"E");

        Double_t QfQbMeanData = QfQbData->GetMean();
        Double_t QfMeanData = QfData->GetMean();
        Double_t QbMeanData = QbData->GetMean();
        Double_t QfQfMeanData = QfQfData->GetMean();

        Double_t bData=(QfQbMeanData-QfMeanData*QbMeanData)/(QfQfMeanData-QfMeanData*QfMeanData);
        HistBData->SetBinContent(i, 14, bData);


        TList *ListMC1 = FileMC->Get(ListNames[0]);
        THnSparseF* SparseMC= (THnSparseF*) ListMC1->FindObject(Form("fPhi2piEta%d",i));
        TH1D *QfQbMCESD = SparseMC->Projection(8,"E");
        TH1D *QfMCESD = SparseMC->Projection(7,"E");
        TH1D *QbMCESD = SparseMC->Projection(6,"E");
        TH1D *QfQfMCESD = SparseMC->Projection(9,"E");

        Double_t QfQbMeanMCESD = QfQbMCESD->GetMean();
        Double_t QfMeanMCESD = QfMCESD->GetMean();
        Double_t QbMeanMCESD = QbMCESD->GetMean();
        Double_t QfQfMeanMCESD = QfQfMCESD->GetMean();

        Double_t bMCESD=(QfQbMeanMCESD-QfMeanMCESD*QbMeanMCESD)/(QfQfMeanMCESD-QfMeanMCESD*QfMeanMCESD);
        HistBMC->SetBinContent(i, 14, bMCESD);


        TList *ListMC2 = FileMC->Get("MCtrue");
        THnSparseF* SparseMCtrue = (THnSparseF*) ListMC2->FindObject(Form("fPhi2piEta%d",i));
        TH1D *QfQbMCtrue = SparseMCtrue->Projection(8,"E");
        TH1D *QfMCtrue = SparseMCtrue->Projection(7,"E");
        TH1D *QbMCtrue = SparseMCtrue->Projection(6,"E");
        TH1D *QfQfMCtrue = SparseMCtrue->Projection(9,"E");

        Double_t QfQbMeanMCtrue = QfQbMCtrue->GetMean();
        Double_t QfMeanMCtrue = QfMCtrue->GetMean();
        Double_t QbMeanMCtrue = QbMCtrue->GetMean();
        Double_t QfQfMeanMCtrue = QfQfMCtrue->GetMean();

        Double_t b_corr=((QfQbMeanData*QfQbMeanMCtrue/QfQbMeanMCESD)-(QfMeanData*QfMeanMCtrue/QfMeanMCESD)*(QbMeanData*QbMeanMCtrue/QbMeanMCESD))/((QfQfMeanData*QfQfMeanMCtrue/QfQfMeanMCESD)-(QfMeanData*QfMeanMCtrue/QfMeanMCESD)*(QfMeanData*QfMeanMCtrue/QfMeanMCESD));

        HistB->SetBinContent(i, 14, b_corr);
    }


    //FB-BF symmetry_____________________________________________________________________________________
    for (Int_t i=1;i<=16;i++)
    {

        TList *ListData2 = FileData->Get(ListNames[0]);
        THnSparseF* SparseData= (THnSparseF*) ListData2->FindObject(Form("fPhi2piEta%d",i));
        TH2D* CorrPlotData=SparseData->Projection(7,6,"");
        TH1D* CorrFuncData=(TH1D*)CorrPlotData->ProfileX("");
        CorrFuncData->Fit("pol1","");
        Double_t bData=CorrFuncData->GetFunction("pol1")->GetParameter(1);
        HistBData->SetBinContent(i,15, bData);


        TList *ListMC3 = FileMC->Get(ListNames[0]);
        THnSparseF* SparseMC= (THnSparseF*) ListMC3->FindObject(Form("fPhi2piEta%d",i));
        TH2D* CorrPlotMC=SparseMC->Projection(7,6,"");
        TH1D* CorrFuncMC=(TH1D*)CorrPlotMC->ProfileX("");
        CorrFuncMC->Fit("pol1","");
        Double_t bMC=CorrFuncMC->GetFunction("pol1")->GetParameter(1);
        HistBMC->SetBinContent(i,15, bMC);


        TList *ListMC4 = FileMC->Get("MCtrue");
        THnSparseF* SparseMCtrue = (THnSparseF*) ListMC4->FindObject(Form("fPhi2piEta%d",i));
        TH2D* CorrPlotMCtrue=SparseMCtrue->Projection(7,6,"");
        TH1D* CorrFuncMCtrue=(TH1D*)CorrPlotMCtrue->ProfileX("");
        CorrFuncMCtrue->Fit("pol1","");
        Double_t bMCtrue=CorrFuncMCtrue->GetFunction("pol1")->GetParameter(1);
        if (bMC){ Double_t b_corr = bData*bMCtrue/bMC;} else {cout<<"ERROR line~175";return;}

        HistB->SetBinContent(i, 15, b_corr);
    }


    //Phojet corrections______________________________________________________________________________
    for (Int_t i=1;i<=16;i++)
    {

        Double_t bData = HistBData->GetBinContent(i,1);
        Double_t bMCtruePhojet=HistBMCtruePhojet->GetBinContent(i);

        TList *ListMC5 = FileMCPhojet->Get(ListNames[0]);
        THnSparseF* SparseMCPhojet= (THnSparseF*) ListMC5->FindObject(Form("fPhi2piEta%d",i));
        TH2D* CorrPlotMCPhojet=SparseMCPhojet->Projection(6,7,"");
        TH1D* CorrFuncMCPhojet=(TH1D*)CorrPlotMCPhojet->ProfileX("");
        CorrFuncMCPhojet->Fit("pol1","");
        Double_t bMCPhojet=CorrFuncMCPhojet->GetFunction("pol1")->GetParameter(1);
        HistBMC->SetBinContent(i,16, bMCPhojet);

        if (bMCPhojet){ Double_t b_corrPhojet = bData*bMCtruePhojet/bMCPhojet;} else {cout<<"ERROR line~215";return;}

        HistB->SetBinContent(i, 16, b_corrPhojet);
    }


    TFile outputfile("Hists.root","RECREATE");
    HistBData->Write();
    HistBMC->Write();
    HistBMCtrue->Write();
    HistBMCtruePhojet->Write();
    HistB->Write();
    outputfile.Close();

 }

